################################################################################
### Main Replication Code                                                    ###
### Title: Do External Threats Increase Bipartisanship in the United States? ###
###        An Experimental Test in the Shadow of China's Rise                ###
### Authors: Eddy S. F. Yeung, Weifang Xu                                    ###
### Version: September 24, 2024                                              ###
################################################################################

### Set-up ----
## Clean the working environment and set up the working directory
rm(list = ls())
# setwd()  # set your working directory here, which should also contain the survey data ("survey_data.csv")

## Load the required packages and functions
library(tidyverse)  # version 2.0.0
library(estimatr)   # version 1.0.0
library(cowplot)    # version 1.1.1
library(texreg)     # version 1.37.5
library(multcomp)   # version 1.4.17
library(interflex)  # version 1.2.6
rescale <- function(x, min, max) {(x - min) / (max - min)}

## Import the dataset
df <- read.csv("survey_data.csv")

### Experimental conditions ----
## Indicators of treatment status
df <- df %>% mutate(exp_condition = case_when(
  exp_3 == 1 ~ "No Threat Prime",
  exp_3 == 2 ~ "Biden Threat Prime",
  exp_3 == 3 ~ "Trump Threat Prime",
  exp_3 == 4 ~ "Nonpartisan Threat Prime"
))
df$exp_condition <- factor(df$exp_condition,
                           levels = c("No Threat Prime", "Biden Threat Prime", 
                                      "Trump Threat Prime", "Nonpartisan Threat Prime"))
table(df$exp_condition)

### Covariates ----
## Age
df$age <- factor(df$age, levels = 1:6,
                 labels = c("18-29", "30-39", "40-49", "50-59", "60-69", "70+"))
sink("demographics.txt")
table(df$age) %>% prop.table() * 100

## Sex
df$sex <- factor(df$sex, levels = 1:2,
                 labels = c("Male", "Female"))
table(df$sex) %>% prop.table() * 100

## Ethnicity
df <- df %>% mutate(eth = case_when(
  race == 1 ~ "White",
  race == 2 ~ "Black",
  race == 4 ~ "Hispanic",
  race == 3 | race == 5 | race == 6 ~ "Other"
))
df$eth <- factor(df$eth, 
                 levels = c("White", "Black", "Hispanic", "Other"))
df$white <- ifelse(df$race == 1, 1, 0)
table(df$eth) %>% prop.table() * 100
sink() # Table S1 is "demographics.tex", which is edited from "demographics.txt" exported here

## Education
df <- df %>% mutate(educ = case_when(
  educ == 1 ~ "No HS",
  educ == 2 ~ "HS",
  educ == 3 ~ "Some college",
  educ == 4 ~ "4-Year College",
  educ == 5 ~ "Post-grad"
))
df$educ <- factor(df$educ, 
                  levels = c("No HS", "HS", "Some college", "4-Year College", "Post-grad"))

## Income
df <- df %>% mutate(income = case_when(
  income == 1 ~ "Less than 30K",
  income == 2 ~ "30K to 70K",
  income == 3 ~ "70K to 100K",
  income == 4 ~ "100K to 200K",
  income == 5 ~ "More than 200K"
))
df$income <- factor(df$income,
                    levels = c("Less than 30K", "30K to 70K", "70K to 100K", "100K to 200K", "More than 200K"))

## Partisanship (-3 = strong Democrat; -3 = strong Republican)
df <- df %>% mutate(pid7 = case_when(
  pid_1 == 2 & pid_2d == 1 ~ -3,
  pid_1 == 2 & pid_2d == 2 ~ -2,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 2 ~ -1,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 4 ~ 0,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 1 ~ 1,
  pid_1 == 1 & pid_2r == 2 ~ 2,
  pid_1 == 1 & pid_2r == 1 ~ 3
))

df <- df %>% mutate(pid3 = case_when(
  pid7 ==  0 ~ "Independent",
  pid7 == -3 | pid7 == -2 | pid7 == -1 ~ "Democrat",
  pid7 ==  3 | pid7 ==  2 | pid7 ==  1 ~ "Republican"
))
df$pid3 <- factor(df$pid3,
                  levels = c("Independent", "Democrat", "Republican"))

## Self-reported ideology (-3 = extremely liberal; 3 = extremely conservative)
df$ideo <- df$ideo - 4

## Nationalism (0 = lowest; 4 = highest)
df$nationalism <- 5 - df$nationalism

## Patriotism (0 = lowest; 3 = highest)
df$patriotism <- 4 - df$patriotism

## National identity (0 = weakest; 4 = strongest)
df$nat_id <- 5 - df$nat_id

## Cooperative internationalism (0 = min; 1 = max)
df$coop_int <- (df$coop_int_1 + df$coop_int_2 + df$coop_int_3 + df$coop_int_4) / 4
df$coop_int <- rescale(df$coop_int, 1, 5)

### Dependent variables ----
## Support for expanding military spending
df$mil_spending <- df$foreign_policy_1 - 1

## Support for restricting scientific exchange
df$sci_exchange <- df$foreign_policy_2 - 1

## Support for reducing bilateral trade
df$trade_reduce <- df$foreign_policy_3 - 1

## Support for offering industry subsidy
df$industry_sub <- df$foreign_policy_4 - 1

## Feeling thermometer (0 = least polarized; 1 = most polarized)
df <- df %>% mutate(affpol_FT = case_when(
  pid3 == "Democrat" ~ affpol_FT_1 - affpol_FT_2,
  pid3 == "Republican" ~ affpol_FT_2 - affpol_FT_1
))
df$affpol_FT <- rescale(df$affpol_FT, -100, 100)

## Trait ratings (0 = most positive; 1 = most negative)
df$affpol_trait_1 <- 6 - df$affpol_trait_1
df$affpol_trait_2 <- 6 - df$affpol_trait_2
df$affpol_trait_3 <- 6 - df$affpol_trait_3
df$affpol_trait_4 <- 6 - df$affpol_trait_4
df$affpol_trait_5 <- 6 - df$affpol_trait_5
df$affpol_trait <- (df$affpol_trait_1 + df$affpol_trait_2 + df$affpol_trait_3 + 
                      df$affpol_trait_4 + df$affpol_trait_5 + df$affpol_trait_6 + 
                      df$affpol_trait_7 + df$affpol_trait_8) / 8
df$affpol_trait <- rescale(df$affpol_trait, 1, 5)

## Trust ratings (0 = most trust; 1 = least trust)
df <- df %>% mutate(affpol_trust = case_when(
  pid3 == "Democrat" ~ 6 - affpol_trust_gop,
  pid3 == "Republican" ~ 6 - affpol_trust_dem
))
df$affpol_trust <- rescale(df$affpol_trust, 1, 5)

## Social distance (0 = least distance; 1 = most distance)
df$affpol_dist_1 <- 5 - df$affpol_dist_1
df$affpol_dist_2 <- 5 - df$affpol_dist_2
df$affpol_dist <- (df$affpol_dist_1 + df$affpol_dist_2 + df$affpol_dist_3) / 3
df$affpol_dist <- rescale(df$affpol_dist, 1, 4)

### Statistical analyses for ATEs on foreign policy preferences ----
## ATE on support for expanding military spending
reg_mil_spending_ATE <- 
  lm_robust(mil_spending ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on restricting scientific exchange
reg_sci_exchange_ATE <- 
  lm_robust(sci_exchange ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on reducing bilateral trade
reg_trade_reduce_ATE <- 
  lm_robust(trade_reduce ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on offering industry subsidy
reg_industry_sub_ATE <- 
  lm_robust(industry_sub ~ exp_condition, data = df, subset = pid3 != "Independent")

### Figure 1: Average Treatment Effects on Preferences for US Foreign Policy ----
ATE_policy <- 
  bind_rows(list(tidy(reg_mil_spending_ATE), tidy(reg_sci_exchange_ATE),
                 tidy(reg_trade_reduce_ATE), tidy(reg_industry_sub_ATE)))
ATE_policy <- subset(ATE_policy, term != "(Intercept)")
ATE_policy <- ATE_policy %>% mutate(treatment = case_when(
  term == "exp_conditionBiden Threat Prime" ~ "Biden Threat Prime",
  term == "exp_conditionTrump Threat Prime" ~ "Trump Threat Prime",
  term == "exp_conditionNonpartisan Threat Prime" ~ "Nonpartisan Threat Prime"
))
ATE_policy$treatment <- 
  factor(ATE_policy$treatment,
         levels = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"))
ATE_policy$outcome <- 
  factor(ATE_policy$outcome,
         levels = c("mil_spending", "sci_exchange", "trade_reduce", "industry_sub"),
         labels = c("Expanding Military Spending", "Restricting Scientific Exchange",
                    "Reducing U.S.-China Trade", "Offering Industry Subsidy"))
ggplot(data = ATE_policy, 
       aes(x = treatment, y = estimate, color = outcome, shape = outcome)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("Foreign Policy", values = c("grey0", "grey25", "grey50", "grey75")) +
  scale_shape_manual("Foreign Policy", values = c(19, 15, 17, 4)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("Average Treatment Effect") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = c(.6985, .896),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.2, 0.6)) + 
  guides(color = guide_legend(nrow = 2, byrow = T))
ggsave(file = "ATE_policy.pdf", width = 8, height = 4)

### Figures S2-S5: Mean and Distribution of Preferences for Foreign Policy Preferences ----
## Support for expanding military spending
summary_temp_dem <- df %>% 
  subset(pid3 == "Democrat") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(mil_spending ~ 1, data = .))) %>% 
  mutate(mil_spending = estimate, partisanship = "Democrat")
summary_temp_gop <- df %>% 
  subset(pid3 == "Republican") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(mil_spending ~ 1, data = .))) %>% 
  mutate(mil_spending = estimate, partisanship = "Republican")
summary_temp <- bind_rows(summary_temp_dem, summary_temp_gop)

p1 <- 
  ggplot(summary_temp, aes(x = exp_condition, y = mil_spending)) + 
  geom_jitter(data = subset(df, pid3 == "Democrat"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#000099", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_jitter(data = subset(df, pid3 == "Republican"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#990000", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_point(aes(color = partisanship, shape = partisanship), 
             size = 2, position = position_dodge(0.15)) +
  scale_color_manual(values = c("#000099", "#990000")) +
  scale_shape_manual(values = c(19, 15)) +
  geom_errorbar(aes(color = partisanship, ymin = conf.low, ymax = conf.high),
                linewidth = 0.5, width = 0, position = position_dodge(.15)) +
  theme_bw() +
  xlab("") +
  ylab("Support for Expanding Military Spending") +
  coord_cartesian(ylim = c(-0.15, 4.15)) +
  theme(text = element_text(color = "black", family = "Times", size = 12),
        axis.text = element_text(color = "black", family = "Times", size = 12),
        legend.justification = c(1, 1), legend.position = c(1, .12),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(color = NA, fill = "transparent"), 
        legend.key = element_rect(color = "transparent", fill = "transparent"),
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank())
ggsave("distribution_policy1.pdf", p1, width = 8, height = 5)

## Support for restricting scientific exchange
summary_temp_dem <- df %>% 
  subset(pid3 == "Democrat") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(sci_exchange ~ 1, data = .))) %>% 
  mutate(sci_exchange = estimate, partisanship = "Democrat")
summary_temp_gop <- df %>% 
  subset(pid3 == "Republican") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(sci_exchange ~ 1, data = .))) %>% 
  mutate(sci_exchange = estimate, partisanship = "Republican")
summary_temp <- bind_rows(summary_temp_dem, summary_temp_gop)

p2 <- 
  ggplot(summary_temp, aes(x = exp_condition, y = sci_exchange)) + 
  geom_jitter(data = subset(df, pid3 == "Democrat"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#000099", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_jitter(data = subset(df, pid3 == "Republican"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#990000", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_point(aes(color = partisanship, shape = partisanship), 
             size = 2, position = position_dodge(0.15)) +
  scale_color_manual(values = c("#000099", "#990000")) +
  scale_shape_manual(values = c(19, 15)) +
  geom_errorbar(aes(color = partisanship, ymin = conf.low, ymax = conf.high),
                linewidth = 0.5, width = 0, position = position_dodge(.15)) +
  theme_bw() +
  xlab("") +
  ylab("Support for Restricting Scientific Exchange") +
  coord_cartesian(ylim = c(-0.15, 4.15)) +
  theme(text = element_text(color = "black", family = "Times", size = 12),
        axis.text = element_text(color = "black", family = "Times", size = 12),
        legend.justification = c(1, 1), legend.position = c(1, .12),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(color = NA, fill = "transparent"), 
        legend.key = element_rect(color = "transparent", fill = "transparent"),
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank())
ggsave("distribution_policy2.pdf", p2, width = 8, height = 5)

## Support for reducing bilateral trade
summary_temp_dem <- df %>% 
  subset(pid3 == "Democrat") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(trade_reduce ~ 1, data = .))) %>% 
  mutate(trade_reduce = estimate, partisanship = "Democrat")
summary_temp_gop <- df %>% 
  subset(pid3 == "Republican") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(trade_reduce ~ 1, data = .))) %>% 
  mutate(trade_reduce = estimate, partisanship = "Republican")
summary_temp <- bind_rows(summary_temp_dem, summary_temp_gop)

p3 <- 
  ggplot(summary_temp, aes(x = exp_condition, y = trade_reduce)) + 
  geom_jitter(data = subset(df, pid3 == "Democrat"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#000099", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_jitter(data = subset(df, pid3 == "Republican"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#990000", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_point(aes(color = partisanship, shape = partisanship), 
             size = 2, position = position_dodge(0.15)) +
  scale_color_manual(values = c("#000099", "#990000")) +
  scale_shape_manual(values = c(19, 15)) +
  geom_errorbar(aes(color = partisanship, ymin = conf.low, ymax = conf.high),
                linewidth = 0.5, width = 0, position = position_dodge(.15)) +
  theme_bw() +
  xlab("") +
  ylab("Support for Reducing Bilateral Trade") +
  coord_cartesian(ylim = c(-0.15, 4.15)) +
  theme(text = element_text(color = "black", family = "Times", size = 12),
        axis.text = element_text(color = "black", family = "Times", size = 12),
        legend.justification = c(1, 1), legend.position = c(1, .12),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(color = NA, fill = "transparent"), 
        legend.key = element_rect(color = "transparent", fill = "transparent"),
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank())
ggsave("distribution_policy3.pdf", p3, width = 8, height = 5)

## Support for offering industry subsidy
summary_temp_dem <- df %>% 
  subset(pid3 == "Democrat") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(industry_sub ~ 1, data = .))) %>% 
  mutate(industry_sub = estimate, partisanship = "Democrat")
summary_temp_gop <- df %>% 
  subset(pid3 == "Republican") %>% 
  group_by(exp_condition) %>%
  do(tidy(lm_robust(industry_sub ~ 1, data = .))) %>% 
  mutate(industry_sub = estimate, partisanship = "Republican")
summary_temp <- bind_rows(summary_temp_dem, summary_temp_gop)

p4 <- 
  ggplot(summary_temp, aes(x = exp_condition, y = industry_sub)) + 
  geom_jitter(data = subset(df, pid3 == "Democrat"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#000099", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_jitter(data = subset(df, pid3 == "Republican"), 
              size = 1, alpha = 0.05, na.rm = T, color = "#990000", 
              position = position_jitter(width = 0.1, height = 0.2)) +
  geom_point(aes(color = partisanship, shape = partisanship), 
             size = 2, position = position_dodge(0.15)) +
  scale_color_manual(values = c("#000099", "#990000")) +
  scale_shape_manual(values = c(19, 15)) +
  geom_errorbar(aes(color = partisanship, ymin = conf.low, ymax = conf.high),
                linewidth = 0.5, width = 0, position = position_dodge(.15)) +
  theme_bw() +
  xlab("") +
  ylab("Support for Offering Industry Subsidy") +
  coord_cartesian(ylim = c(-0.15, 4.15)) +
  theme(text = element_text(color = "black", family = "Times", size = 12),
        axis.text = element_text(color = "black", family = "Times", size = 12),
        legend.justification = c(1, 1), legend.position = c(1, .12),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(color = NA, fill = "transparent"), 
        legend.key = element_rect(color = "transparent", fill = "transparent"),
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank())
ggsave("distribution_policy4.pdf", p4, width = 8, height = 5)

### Statistical analyses for ATEs on affective polarization ----
## ATE on feeling thermometer
reg_affpol_FT_ATE <- 
  lm_robust(affpol_FT ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on trait ratings
reg_affpol_trait_ATE <- 
  lm_robust(affpol_trait ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on trust ratings
reg_affpol_trust_ATE <- 
  lm_robust(affpol_trust ~ exp_condition, data = df, subset = pid3 != "Independent")

## ATE on social distance
reg_affpol_dist_ATE <- 
  lm_robust(affpol_dist ~ exp_condition, data = df, subset = pid3 != "Independent")

### Figure 3: Average Treatment Effects on Affective Polarization ----
ATE_affpol <- 
  bind_rows(list(tidy(reg_affpol_FT_ATE), tidy(reg_affpol_trait_ATE), 
                 tidy(reg_affpol_trust_ATE), tidy(reg_affpol_dist_ATE)))
ATE_affpol <- subset(ATE_affpol, term != "(Intercept)")
ATE_affpol <- ATE_affpol %>% mutate(treatment = case_when(
  term == "exp_conditionBiden Threat Prime" ~ "Biden Threat Prime",
  term == "exp_conditionTrump Threat Prime" ~ "Trump Threat Prime",
  term == "exp_conditionNonpartisan Threat Prime" ~ "Nonpartisan Threat Prime"
))
ATE_affpol$treatment <- 
  factor(ATE_affpol$treatment,
         levels = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"))
ATE_affpol$outcome <- 
  factor(ATE_affpol$outcome,
         levels = c("affpol_FT", "affpol_trait", "affpol_trust", "affpol_dist"),
         labels = c("Feeling Thermometer", "Outparty Trait Ratings",
                    "Outparty Trust Ratings", "Outparty Social Distance"))
ggplot(data = ATE_affpol, 
       aes(x = treatment, y = estimate, color = outcome, shape = outcome)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("Affective Polarization", values = c("grey0", "grey25", "grey50", "grey75")) +
  scale_shape_manual("Affective Polarization", values = c(19, 15, 17, 4)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("Average Treatment Effect") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = c(.748, .896),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.15, 0.15)) + 
  guides(color = guide_legend(nrow = 2, byrow = T))
ggsave(file = "ATE_affpol.pdf", width = 8, height = 4)

### Table S3: OLS Regression Corresponding to Equation (1) ----
df$dem <- ifelse(df$pid3 == "Democrat", 1, 0)

## Support for expanding military spending
mod1.1 <- 
  lm_robust(mil_spending ~ exp_condition * dem, 
            data = df, subset = pid3 != "Independent")
mod1.2 <- 
  lm_robust(mil_spending ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## Support for restricting scientific exchange
mod2.1 <- 
  lm_robust(sci_exchange ~ exp_condition * dem, 
            data = df, subset = pid3 != "Independent")
mod2.2 <- 
  lm_robust(sci_exchange ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## Support for reducing bilateral trade
mod3.1 <- 
  lm_robust(trade_reduce ~ exp_condition * dem, 
            data = df, subset = pid3 != "Independent")
mod3.2 <- 
  lm_robust(trade_reduce ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## Support for offering industry subsidy
mod4.1 <- 
  lm_robust(industry_sub ~ exp_condition * dem, 
            data = df, subset = pid3 != "Independent")
mod4.2 <- 
  lm_robust(industry_sub ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## Make a regression table
sink("table_regression.txt")
texreg(list(mod1.1, mod1.2, mod2.1, mod2.2, mod3.1, mod3.2, mod4.1, mod4.2),
       stars = c(0.001, 0.01, 0.05), 
       include.ci = F,
       custom.header = list("Dependent Variable: 5-Point Support for the Policy" = 1:8),
       custom.model.names = c("Military Spending", "Military Spending", 
                              "Scientific Exchange", "Scientific Exchange",
                              "Trade Reduction", "Trade Reduction", 
                              "Industry Subsidy", "Industry Subsidy"),
       custom.coef.names = 
         c("Constant", "Biden Threat Prime", "Trump Threat Prime", 
           "Nonpartisan Threat Prime", "Democrat", 
           "Biden Threat Prime × Democrat", "Trump Threat Prime × Democrat",
           "Nonpartisan Threat Prime × Democrat", "Aged 30-39",
           "Aged 40-49", "Aged 50-59", "Aged 60-69", "Aged 70+",
           "Female", "White", "High School", "Some College", "4-Year College", 
           "Postgraduate Degree", "Income 30K-70K", "Income 70K-100K", 
           "Income 100K-200K", "Income 200K+", "Self-Reported Ideology (7-point)",
           "Nationalism (5-point)", "Patriotism (4-point)", 
           "National Identity (5-point)", "Cooperative Internationalism"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.",
       caption = "OLS Regression Corresponding to Equation 1",
       fontsize = "small"
)
sink() # Table S3 is "table_regression.tex", which is edited from "table_regression.txt" exported here

### Figure 2: CATEs and Difference-in-CATEs for Democrats' and Republicans' US Foreign Policy Preferences ----
## Combine the estimates to one data frame and keep relevant estimates only
CATE_summary_df <- 
  bind_rows(tidy(mod1.2), tidy(mod2.2), tidy(mod3.2), tidy(mod4.2)) %>% 
  filter(term == "exp_conditionBiden Threat Prime" | 
           term == "exp_conditionTrump Threat Prime" | 
           term == "exp_conditionNonpartisan Threat Prime" | 
           term == "exp_conditionBiden Threat Prime:dem" | 
           term == "exp_conditionTrump Threat Prime:dem" | 
           term == "exp_conditionNonpartisan Threat Prime:dem")

## Rename the terms
CATE_summary_df$term <- 
  factor(CATE_summary_df$term, 
         levels = c("exp_conditionBiden Threat Prime", 
                    "exp_conditionTrump Threat Prime", 
                    "exp_conditionNonpartisan Threat Prime",
                    "exp_conditionBiden Threat Prime:dem",
                    "exp_conditionTrump Threat Prime:dem",
                    "exp_conditionNonpartisan Threat Prime:dem"),
         labels = c(rep("CATE on Republicans", 3),
                    rep("Difference in CATEs (CATE on Dems - CATE on Reps)", 3)))

## Indicate treatment conditions for the estimates
CATE_summary_df <- CATE_summary_df %>% mutate(treatment = rep(1:3, 8))
CATE_summary_df$treatment <- 
  factor(CATE_summary_df$treatment, levels = 1:3,
         labels = c("Biden\nThreat Prime", "Trump\nThreat Prime", "Nonpartisan\nThreat Prime"))

## Estimate the CATEs on Democrats and merge them to the data frame
# Military spending (Biden treatment)
temp_estimate_1 <- glht(mod1.2, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Scientific exchange (Biden treatment)
temp_estimate_2 <- glht(mod2.2, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Bilateral trade (Biden treatment)
temp_estimate_3 <- glht(mod3.2, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Industry subsidy (Biden treatment)
temp_estimate_4 <- glht(mod4.2, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Merge all estimates for Biden treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)

# Military spending (Trump treatment)
temp_estimate_1 <- glht(mod1.2, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Scientific exchange (Trump treatment)
temp_estimate_2 <- glht(mod2.2, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Bilateral trade (Trump treatment)
temp_estimate_3 <- glht(mod3.2, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Industry subsidy (Trump treatment)
temp_estimate_4 <- glht(mod4.2, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Merge all estimates for Trump treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)

# Military spending (Nonpartisan treatment)
temp_estimate_1 <- glht(mod1.2, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Scientific exchange (Nonpartisan treatment)
temp_estimate_2 <- glht(mod2.2, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Bilateral trade (Nonpartisan treatment)
temp_estimate_3 <- glht(mod3.2, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Industry subsidy (Nonpartisan treatment)
temp_estimate_4 <- glht(mod4.2, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

# Merge all estimates for Nonpartisan treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)
CATE_summary_df$term <- relevel(CATE_summary_df$term, "CATE on Democrats")
CATE_summary_df$outcome <- 
  factor(CATE_summary_df$outcome, 
         levels = c("mil_spending", "sci_exchange", "trade_reduce", "industry_sub"),
         labels = c("Expanding Military Spending", "Restricting Scientific Exchange",
                    "Reducing U.S.-China Trade", "Offering Industry Subsidy"))

## Plot the estimates
ggplot(data = CATE_summary_df, 
       aes(x = treatment, y = estimate, color = term, shape = term)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("", values = c("grey0", "grey40", "grey80")) +
  scale_shape_manual("", values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("Estimated Treatment Effect") +
  facet_wrap(~outcome, nrow = 2) +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.5, 0.5)) + 
  guides(color = guide_legend(nrow = 1, byrow = T))
ggsave("CATE_policy.pdf", width = 8, height = 6)

### Exploratory analyses for Independents ----
## Sample size per group
table(df$pid3, df$exp_condition)

## Policy preferences
# CATE on Independents' support for expanding military spending
reg_mil_spending_CATE_ind <- 
  lm_robust(mil_spending ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' support for restricting scientific exchange
reg_sci_exchange_CATE_ind <-
  lm_robust(sci_exchange ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' support for reducing bilateral trade
reg_trade_reduce_CATE_ind <- 
  lm_robust(trade_reduce ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' support for offering industry subsidy
reg_industry_sub_CATE_ind <- 
  lm_robust(industry_sub ~ exp_condition, data = df, subset = pid3 == "Independent")

# Figure S14: CATEs for Independents' Foreign Policy Preferences ----
CATE_policy_ind <- 
  bind_rows(list(tidy(reg_mil_spending_CATE_ind), tidy(reg_sci_exchange_CATE_ind),
                 tidy(reg_trade_reduce_CATE_ind), tidy(reg_industry_sub_CATE_ind)))
CATE_policy_ind <- subset(CATE_policy_ind, term != "(Intercept)")
CATE_policy_ind <- CATE_policy_ind %>% mutate(treatment = case_when(
  term == "exp_conditionBiden Threat Prime" ~ "Biden Threat Prime",
  term == "exp_conditionTrump Threat Prime" ~ "Trump Threat Prime",
  term == "exp_conditionNonpartisan Threat Prime" ~ "Nonpartisan Threat Prime"
))
CATE_policy_ind$treatment <- 
  factor(CATE_policy_ind$treatment,
         levels = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"))
CATE_policy_ind$outcome <- 
  factor(CATE_policy_ind$outcome,
         levels = c("mil_spending", "sci_exchange", "trade_reduce", "industry_sub"),
         labels = c("Expanding Military Spending", "Restricting Scientific Exchange",
                    "Reducing U.S.-China Trade", "Offering Industry Subsidy"))
ggplot(data = CATE_policy_ind, 
       aes(x = treatment, y = estimate, color = outcome, shape = outcome)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("Foreign Policy", values = c("grey0", "grey25", "grey50", "grey75")) +
  scale_shape_manual("Foreign Policy", values = c(19, 15, 17, 4)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Independents") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = c(.698, .896),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.5, 1)) + 
  guides(color = guide_legend(nrow = 2, byrow = T))
ggsave(file = "CATE_policy_ind.pdf", width = 8, height = 4)

## Party affect
# CATE on Independents' feelings toward Democrats
df$affect_dem <- rescale(df$affpol_FT_1, 0, 100)
reg_affect_dem_CATE_ind <- 
  lm_robust(affect_dem ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' feelings toward Republicans
df$affect_gop <- rescale(df$affpol_FT_2, 0, 100)
reg_affect_gop_CATE_ind <- 
  lm_robust(affect_gop ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' trust in the Democratic Party
df$trust_dem <- rescale(df$affpol_trust_dem, 1, 5)
reg_trust_dem_CATE_ind <- 
  lm_robust(trust_dem ~ exp_condition, data = df, subset = pid3 == "Independent")

# CATE on Independents' trust in the Republican Party
df$trust_gop <- rescale(df$affpol_trust_gop, 1, 5)
reg_trust_gop_CATE_ind <- 
  lm_robust(trust_gop ~ exp_condition, data = df, subset = pid3 == "Independent")

# Figure S15: CATEs for Independents’ toward Democrats and Republicans ----
CATE_affect_ind <- 
  bind_rows(list(tidy(reg_affect_dem_CATE_ind), tidy(reg_affect_gop_CATE_ind), 
                 tidy(reg_trust_dem_CATE_ind), tidy(reg_trust_gop_CATE_ind)))
CATE_affect_ind <- subset(CATE_affect_ind, term != "(Intercept)")
CATE_affect_ind <- CATE_affect_ind %>% mutate(treatment = case_when(
  term == "exp_conditionBiden Threat Prime" ~ "Biden Threat Prime",
  term == "exp_conditionTrump Threat Prime" ~ "Trump Threat Prime",
  term == "exp_conditionNonpartisan Threat Prime" ~ "Nonpartisan Threat Prime"
))
CATE_affect_ind$treatment <- 
  factor(CATE_affect_ind$treatment,
         levels = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"))
CATE_affect_ind$outcome <- 
  factor(CATE_affect_ind$outcome,
         levels = c("affect_dem", "affect_gop", "trust_dem", "trust_gop"),
         labels = c("Feelings toward Democrats", "Feelings toward Republicans",
                    "Trust in Democratic Party", "Trust in Republican Party"))
ggplot(data = CATE_affect_ind, 
       aes(x = treatment, y = estimate, color = outcome, shape = outcome)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("Party Affect", values = c("grey0", "grey25", "grey50", "grey75")) +
  scale_shape_manual("Party Affect", values = c(19, 15, 17, 4)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Independents") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = c(.717, .896),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.15, 0.15)) + 
  guides(color = guide_legend(nrow = 2, byrow = T))
ggsave(file = "CATE_affect_ind.pdf", width = 8, height = 4)

### Figures S6-S13: Heterogeneous Treatment Effects on Foreign Policy Preferences ----
## By nationalism
# Expanding military spending
HTE_nat_mil <- 
  interflex(Y = "mil_spending", D = "exp_condition", X = "nationalism",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "patriotism", "nat_id", "coop_int"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_nat_mil, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Nationalism (0 = least nationalistic; 4 = most nationalistic)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_nat_mil.pdf", width = 8, height = 4)

# Restricting scientific exchange
HTE_nat_sci <- 
  interflex(Y = "sci_exchange", D = "exp_condition", X = "nationalism",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "patriotism", "nat_id", "coop_int"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_nat_sci, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Nationalism (0 = least nationalistic; 4 = most nationalistic)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_nat_sci.pdf", width = 8, height = 4)

# Reducing bilateral trade
HTE_nat_trade <- 
  interflex(Y = "trade_reduce", D = "exp_condition", X = "nationalism",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "patriotism", "nat_id", "coop_int"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_nat_trade, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Nationalism (0 = least nationalistic; 4 = most nationalistic)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_nat_trade.pdf", width = 8, height = 4)

# Offering industry subsidy
HTE_nat_sub <- 
  interflex(Y = "industry_sub", D = "exp_condition", X = "nationalism",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "patriotism", "nat_id", "coop_int"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_nat_sub, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Nationalism (0 = least nationalistic; 4 = most nationalistic)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_nat_sub.pdf", width = 8, height = 4)

## By cooperative internationalism
# Expanding military spending
HTE_coop_mil <- 
  interflex(Y = "mil_spending", D = "exp_condition", X = "coop_int",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "nationalism", "patriotism", "nat_id"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_coop_mil, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Cooperative Internationalism (0 = least cooperative; 1 = most cooperative)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_coop_mil.pdf", width = 8, height = 4)

# Restricting scientific exchange
HTE_coop_sci <- 
  interflex(Y = "sci_exchange", D = "exp_condition", X = "coop_int",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "nationalism", "patriotism", "nat_id"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_coop_sci, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Cooperative Internationalism (0 = least cooperative; 1 = most cooperative)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_coop_sci.pdf", width = 8, height = 4)

# Reducing bilateral trade
HTE_coop_trade <- 
  interflex(Y = "trade_reduce", D = "exp_condition", X = "coop_int",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "nationalism", "patriotism", "nat_id"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_coop_trade, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Cooperative Internationalism (0 = least cooperative; 1 = most cooperative)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_coop_trade.pdf", width = 8, height = 4)

# Offering industry subsidy
HTE_coop_sub <- 
  interflex(Y = "industry_sub", D = "exp_condition", X = "coop_int",
            Z = c("age", "sex", "white", "educ", "income", "ideo", "pid7",
                  "nationalism", "patriotism", "nat_id"),
            base = "No Threat Prime",
            data = df, na.rm = TRUE,
            estimator = "binning", vcov.type = "robust")
plot(HTE_coop_sub, order = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"),
     subtitles = c("Biden Threat Prime", "Trump Threat Prime", "Nonpartisan Threat Prime"), 
     theme.bw = TRUE, cex.axis = 1.2, cex.lab = 0.8, bin.labs = TRUE,
     xlab = "Cooperative Internationalism (0 = least cooperative; 1 = most cooperative)",
     ylab = "Marginal Treatment Effect", ylim = c(-0.7, 0.7),
     file = "HTE_coop_sub.pdf", width = 8, height = 4)

### Figure S1: CATEs and Differences-in-CATEs for Democrats' and Republicans' US Foreign Policy Preferences from an Alternative Operationalization ----
## Create ternary variables for foreign policy preferences
df <- df %>% mutate(mil_spending_alt = case_when(
  mil_spending == 0 | mil_spending == 1 ~ -1,
  mil_spending == 2 ~ 0,
  mil_spending == 3 | mil_spending == 4 ~ 1
))
table(df$mil_spending_alt)

df <- df %>% mutate(sci_exchange_alt = case_when(
  sci_exchange == 0 | sci_exchange == 1 ~ -1,
  sci_exchange == 2 ~ 0,
  sci_exchange == 3 | sci_exchange == 4 ~ 1
))
table(df$sci_exchange_alt)

df <- df %>% mutate(trade_reduce_alt = case_when(
  trade_reduce == 0 | trade_reduce == 1 ~ -1,
  trade_reduce == 2 ~ 0,
  trade_reduce == 3 | trade_reduce == 4 ~ 1
))
table(df$trade_reduce_alt)

df <- df %>% mutate(industry_sub_alt = case_when(
  industry_sub == 0 | industry_sub == 1 ~ -1,
  industry_sub == 2 ~ 0,
  industry_sub == 3 | industry_sub == 4 ~ 1
))
table(df$industry_sub_alt)

## HTE on support for expanding military spending
mod1 <- 
  lm_robust(mil_spending_alt ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## HTE on support for restricting scientific exchange
mod2 <- 
  lm_robust(sci_exchange_alt ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## HTE on support for reducing bilateral trade
mod3 <- 
  lm_robust(trade_reduce_alt ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## HTE on support for offering industry subsidy
mod4 <- 
  lm_robust(industry_sub_alt ~ exp_condition * dem + age + sex + white + educ +
              income + ideo + nationalism + patriotism + nat_id + coop_int, 
            data = df, subset = pid3 != "Independent")

## Combine the estimates to one data frame and keep relevant estimates only
CATE_summary_df <- 
  bind_rows(tidy(mod1), tidy(mod2), tidy(mod3), tidy(mod4)) %>% 
  filter(term == "exp_conditionBiden Threat Prime" | 
           term == "exp_conditionTrump Threat Prime" | 
           term == "exp_conditionNonpartisan Threat Prime" | 
           term == "exp_conditionBiden Threat Prime:dem" | 
           term == "exp_conditionTrump Threat Prime:dem" | 
           term == "exp_conditionNonpartisan Threat Prime:dem")

## Rename the terms
CATE_summary_df$term <- 
  factor(CATE_summary_df$term, 
         levels = c("exp_conditionBiden Threat Prime", 
                    "exp_conditionTrump Threat Prime", 
                    "exp_conditionNonpartisan Threat Prime",
                    "exp_conditionBiden Threat Prime:dem",
                    "exp_conditionTrump Threat Prime:dem",
                    "exp_conditionNonpartisan Threat Prime:dem"),
         labels = c(rep("CATE on Republicans", 3),
                    rep("Difference in CATEs (CATE on Dems - CATE on Reps)", 3)))

## Indicate treatment conditions for the estimates
CATE_summary_df <- CATE_summary_df %>% mutate(treatment = rep(1:3, 8))
CATE_summary_df$treatment <- 
  factor(CATE_summary_df$treatment, levels = 1:3,
         labels = c("Biden\nThreat Prime", "Trump\nThreat Prime", "Nonpartisan\nThreat Prime"))

## Estimate the CATEs on Democrats and merge them to the data frame
temp_estimate_1 <- glht(mod1, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending_alt", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_2 <- glht(mod2, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange_alt", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_3 <- glht(mod3, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce_alt", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_4 <- glht(mod4, linfct = c("`exp_conditionBiden Threat Prime` + `exp_conditionBiden Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub_alt", "Biden\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

## Merge all estimates for Biden treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)

temp_estimate_1 <- glht(mod1, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending_alt", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_2 <- glht(mod2, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange_alt", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_3 <- glht(mod3, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce_alt", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_4 <- glht(mod4, linfct = c("`exp_conditionTrump Threat Prime` + `exp_conditionTrump Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub_alt", "Trump\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

## Merge all estimates for Trump treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)

temp_estimate_1 <- glht(mod1, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "mil_spending_alt", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_2 <- glht(mod2, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "sci_exchange_alt", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_3 <- glht(mod3, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "trade_reduce_alt", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

temp_estimate_4 <- glht(mod4, linfct = c("`exp_conditionNonpartisan Threat Prime` + `exp_conditionNonpartisan Threat Prime:dem` = 0"))
temp_df_4 <- 
  as.data.frame(matrix(c("CATE on Democrats", 
                         confint(temp_estimate_4)$confint[ , c("Estimate", "lwr", "upr")], 
                         "industry_sub_alt", "Nonpartisan\nThreat Prime"), 
                       nrow = 1))
temp_df_4 <- temp_df_4 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)

## Merge all estimates for Nonpartisan treatment to the main data frame
CATE_dem_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3, temp_df_4)
CATE_dem_only$term <- as.factor(CATE_dem_only$term)
CATE_dem_only$estimate <- as.numeric(CATE_dem_only$estimate)
CATE_dem_only$conf.low <- as.numeric(CATE_dem_only$conf.low)
CATE_dem_only$conf.high <- as.numeric(CATE_dem_only$conf.high)
CATE_dem_only$outcome <- as.factor(CATE_dem_only$outcome)
CATE_dem_only$treatment <- as.factor(CATE_dem_only$treatment)
CATE_summary_df <- bind_rows(CATE_summary_df, CATE_dem_only)
CATE_summary_df$term <- relevel(CATE_summary_df$term, "CATE on Democrats")
CATE_summary_df$outcome <- 
  factor(CATE_summary_df$outcome, 
         levels = c("mil_spending_alt", "sci_exchange_alt", "trade_reduce_alt", "industry_sub_alt"),
         labels = c("Expanding Military Spending", "Restricting Scientific Exchange",
                    "Reducing U.S.-China Trade", "Offering Industry Subsidy"))

## Plot the estimates
ggplot(data = CATE_summary_df, 
       aes(x = treatment, y = estimate, color = term, shape = term)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual("", values = c("grey0", "grey40", "grey80")) +
  scale_shape_manual("", values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("Estimated Treatment Effect") +
  facet_wrap(~outcome, nrow = 2) +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-0.3, 0.3)) + 
  guides(color = guide_legend(nrow = 1, byrow = T))
ggsave("CATE_policy_alt.pdf", width = 8, height = 6)
